using chc 6
present CV LDAs, have testing/training if need be
graphed LDAs are all the points.
I cannot do site contrasts because of our uneven samplign scheme. The MANOVA throws out samples.
fit_full_species_man$pca_summary
## Importance of first k=7 (out of 35) components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.0736 0.7112 0.47335 0.37674 0.33980 0.30966 0.2955
## Proportion of Variance 0.4066 0.1784 0.07903 0.05006 0.04073 0.03382 0.0308
## Cumulative Proportion 0.4066 0.5850 0.66402 0.71409 0.75481 0.78864 0.8194
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 4 0.03980 40.088 28 776.62 < 2.2e-16 ***
## sex 1 0.81535 6.956 7 215.00 1.795e-07 ***
## hostRace:sex 4 0.79391 1.833 28 776.62 0.005633 **
## Residuals 221
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Cingulata Cornivora Mendax pom Zepheria
## Cingulata 25 0 0 0 0
## Cornivora 1 5 0 0 0
## Mendax 0 0 60 5 1
## pom 0 0 1 116 0
## Zepheria 0 0 0 1 16
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 9.610390e-01 9.390180e-01 9.273307e-01 9.820323e-01 5.281385e-01
## AccuracyPValue McnemarPValue
## 1.508463e-49 NaN
Collapsed all 35 chc measurements into a value column with a corresponding chc column detailing what chc was measured. Fitted a model to evaluate the efect of host and sex on chc values. The model accounts for uneven site effects and individual effects.
All Bayesian models were created in Stan computational framework (http://mc-stan.org/) accessed with brms package (Bürkner, 2017). To improve convergence and guard against overfitting, we specified mildly informative conservative priors. We ran 4 chains with 2000 interations. We checked for convergence and found all models mixed well. (Pareto K estimates are all good!)
#1 value ~ hostRace * sex + (1|site) +(1|chc) + (1|ID)
mod = readRDS("fit_chc_model1.rds")
plot(mod)
brms::pp_check(mod, resp ="value") #great
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
bayes_R2(mod, summary = TRUE)
## Estimate Est.Error Q2.5 Q97.5
## R2 0.8790622 0.0009515364 0.8771306 0.8808946
tab_model(mod)
| Â | value | |
|---|---|---|
| Predictors | Estimates | CI (95%) |
| Intercept | -0.25 | -0.49 – -0.00 |
| hostRace: Cornivora | -0.13 | -0.29 – 0.03 |
| hostRace: Mendax | 0.17 | 0.09 – 0.26 |
| hostRace: Zepheria | 0.33 | 0.23 – 0.43 |
| hostRace: pom | 0.11 | 0.04 – 0.17 |
| sex: Male | 0.03 | -0.06 – 0.11 |
| hostRaceCornivora.sexMale | -0.14 | -0.35 – 0.06 |
| hostRaceMendax.sexMale | -0.07 | -0.17 – 0.03 |
| hostRaceZepheria.sexMale | -0.01 | -0.14 – 0.12 |
| hostRacepom.sexMale | -0.14 | -0.23 – -0.05 |
| Random Effects | ||
| σ2 | 0.44 | |
| τ00 | 0.07 | |
| ICC | 0.86 | |
| N site | 8 | |
| N ID | 231 | |
| N chc | 35 | |
| Observations | 8085 | |
| Marginal R2 / Conditional R2 | 0.021 / 0.879 | |
plot_model(mod, type = "pred", terms = c("hostRace", "sex"))
plot_model(mod, bpe = "mean", bpe.style = "dot", prob.inner = 0.4, prob.outer = 0.8)
summary(manova(as.matrix(data[,4:cols]) ~ hostRace * site * sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## hostRace 1 0.82295 3.3194 7 108 0.003071 **
## site 1 0.85799 2.5537 7 108 0.017937 *
## sex 1 0.52503 13.9576 7 108 8.859e-13 ***
## hostRace:site 1 0.69396 6.8041 7 108 1.085e-06 ***
## hostRace:sex 1 0.95188 0.7799 7 108 0.605479
## site:sex 1 0.93727 1.0326 7 108 0.412735
## hostRace:site:sex 1 0.90238 1.6691 7 108 0.124190
## Residuals 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Apple Haw
## Apple 25 12
## Haw 14 71
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.786885246 0.503288443 0.703538181 0.855813669 0.680327869
## AccuracyPValue McnemarPValue
## 0.006229331 0.844519267
## Reference
## Prediction Apple_Female Apple_Male Haw_Female Haw_Male
## Apple_Female 23 2 9 1
## Apple_Male 1 1 1 2
## Haw_Female 8 0 32 3
## Haw_Male 1 3 5 30
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 7.049180e-01 5.725131e-01 6.155825e-01 7.840219e-01 3.852459e-01
## AccuracyPValue McnemarPValue
## 8.606186e-13 9.110307e-01
# no interaction because missing males at MtPleasant
summary(manova(as.matrix(data[,4:cols]) ~ sex + site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.74708 0.67708 5 10 0.6507
## site 1 0.83749 0.38808 5 10 0.8461
## Residuals 14
## Reference
## Prediction Zepheria_Female Zepheria_Male
## Zepheria_Female 5 3
## Zepheria_Male 4 5
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.5882353 0.1793103 0.3292472 0.8155630 0.5294118
## AccuracyPValue McnemarPValue
## 0.4062810 1.0000000
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.57727 6.1024 6 50 7.541e-05 ***
## site 2 0.10716 17.1232 12 100 < 2.2e-16 ***
## sex:site 2 0.74705 1.3082 12 100 0.2258
## Residuals 55
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Mendax_Female Mendax_Male
## Mendax_Female 20 15
## Mendax_Male 13 13
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.54098361 0.07072905 0.40849889 0.66935590 0.54098361
## AccuracyPValue McnemarPValue
## 0.55241652 0.85010674
summary(manova(as.matrix(data[,4:cols]) ~ sex * site, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.78806 1.2774 4 19 0.3136059
## site 1 0.38209 7.6817 4 19 0.0007371 ***
## sex:site 1 0.70133 2.0229 4 19 0.1319094
## Residuals 22
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Reference
## Prediction Cingulata_Female Cingulata_Male
## Cingulata_Female 5 7
## Cingulata_Male 7 7
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.46153846 -0.08333333 0.26587122 0.66629178 0.53846154
## AccuracyPValue McnemarPValue
## 0.83731354 1.00000000
There is only one site
# no site because only sampled at one site.
summary(manova(as.matrix(data[,4:cols]) ~ sex, data = data), test = "Wilks")
## Df Wilks approx F num Df den Df Pr(>F)
## sex 1 0.54542 0.83346 2 2 0.5454
## Residuals 3
Cant show a plot
## Reference
## Prediction Cornivora_Female Cornivora_Male
## Cornivora_Female 1 0
## Cornivora_Male 1 3
## Accuracy Kappa AccuracyLower AccuracyUpper AccuracyNull
## 0.8000000 0.5454545 0.2835821 0.9949492 0.6000000
## AccuracyPValue McnemarPValue
## 0.3369600 1.0000000